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dislodging  forces  is  the  generation  of  stagnation  points  in  the  flow.  In  order 
to  examine  further  the  role  played  by  flow  separation  around  the  blunt  body 
segments,  a complex  velocity  potential  is  developed  to  describe  a stationary 
vortex  pair  located  in  the  wake  region  of  the  flow.  It  remains  to  super-impose 
this  vortex  pair  on  the  unseparated  cross  flow  in  order  to  ascertain  the  drag 
forces  which  contribute  to  limb  dislodgement. 
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SUMMARY 


Problem 


Statistical  data  show  that  wlndblast  forces  Increase  with  aircraft 
speed  to  the  point  where  an  overall  five  or  ten  percent  limb-flail  injury 
rate  rises  to  40  percent  or  more.  This  is  far  from  negligible,  but  control 
of  the  forces  that  produce  excessive  motion  of  the  limbs  of  an  ejection 
seat  occupant  can  only  be  achieved  if  we  increase  our  understanding  of  the 
aerodynamic  loading  to  which  a pilot  is  exposed  during  high-speed  eject'  ns. 
Toward  this  end,  a mathematical  model  is  being  developed  which  is  expec  ed 
to  provide  aerodynamic  data  that  can  be  used  as  input  information  to  the 
Aerospace  Medical  Research  Laboratory's  Articulated  Total  Body  Model.  The 
latter  can  then  be  used  to  assess  the  kinematics  of  limb  motion  under  the 
action  of  specific  aerodynamic  forces.  In  this  report,  the  forces  tending 
to  dislodge  the  limbs  of  an  ejection  seat  occupant  from  one  another,  or 
from  a restraining  surface,  are  calculated  in  the  absence  of  flow  separa- 
tion. The  results  are  then  modified  to  show  how  the  more  realistic  physical 
effects  of  separation  of  flow  around  the  blunt  body  segments  may  be  taken 
into  consideration. 

Approach 

Limb-dislodging  forces  in  the  absence  of  flow  separation  are  computed 
by  Integrating  the  pressure  coefficient  which  was  derived  in  an  earlier 
report  for  the  cross-flow  over  a limb-to-limb  or  limb-to-restraining  surface 
contact  configuration.  Modifications  to  include  the  effects  of  senara- 
tion  of  flow  are  then  introduced  by  describing  a distribution  of  inviscid 
vortices  which  originate  from  the  separation  of  shear  lavers,  and  which 
can  be  superimposed  on  the  unseparated  potential  flow  solutions  obtained 
thus  far. 


Results 


For  a simulated  ejection  taking  place  at  an  altitude  of  10,000  feet, 
the  results  show  that,  based  on  published  data  for  a pilot's  average  grip 
retention  capability,  the  probability  of  his  letting  go  is  100%  if  the 
ejection  occurs  at  Mach  numbers  in  excess  of  0.72.  That  is,  the  aerody- 
namic forces  are  such  that  the  pilot's  musculo-skeletal  system  is  not 
likely  to  withstand  the  tendency  for  dislodgement  from  a restraining  sur- 
face if  he  is  ejecting  at  Mach  numbers  exceeding  0.72.  Moreover,  the 
probability  of  major  flail  injurv  is  around  100%  if  the  ejection  Mach 
number  exceeds  1.25.  One  major  factor  which  contributes  to  these  large 
limb-dislodging  forces  appears  to  be  the  generation  of  stagnation  points 
in  the  flow.  In  order  to  examine  further  the  role  played  by  flow  separa- 
tion around  the  blunt  body  segments,  a complex  velocity  potential  is 
developed  to  describe  a stationary  vortex  pair  located  in  the  wake  region 
of  the  flow.  It  remains  to  superimpose  this  vortex  pair  on  the  unseparated 
cross  flow  in  order  to  ascertain  the  drag  forces  which  contribute  to  limb 
dislodgement. 

Conclusions 


To  the  level  of  approximation  in  this  report,  one  would  have  to  con- 
clude that  the  generation  of  stagnation  points  in  the  flow  produces  forces 
that  can  cause  linb-dislodgement  (with  subsequent  flail  and  possible  serious 
injury).  Moreover,  these  forces  are  sensitive  to  the  angle  at  which  the 
limb  Intercepts  the  flow,  such  that  the  higher  the  angle,  the  greater  the 
tendency  for  dislodgement.  And  finally,  the  forces  increase  rapidly  with 
speed  of  ejection,  which  correlates  well  with  the  finding  that  windblast 
injuries  increase  dramatically  as  a function  of  airspeed. 

Recommendat ions 


It  is  desirable  to  examine  the  time-course  of  the  limb  dislodging  forces 
after  the  onset  of  windblast.  That  is,  the  forces  computed  in  this  report 
exist  at  the  onset  of  ejection  and  are  peak  forces  that  prevail  so  long 
as  the  limb  is  in  contact  with  a restraining  surface.  Once  contact  with 
the  surface  is  broken,  there  follows  a redistribution  of  forces,  and  this 
information  is  desirable  input  to  the  ATB  simulation  which  provides  kine- 
matic data.  The  added  complications  resulting  from  separation  of  the 
flow  and  limb  Interactions  need  also  to  be  further  elucidated  in  an  effort 
to  define  design  criteria  for  safer  ejections.  There  is  interest  in  extend- 
ing the  analysis  to  include  the  head/neck  configuration  in  order  to  shed 
some  light  on  the  factors  that  contribute  to  an  oscillatory  fore-aft  head 
vibration.  And  finally,  the  effects  of  surface  roughness  (which  may  be 
due,  in  part,  to  the  type  of  clothing  being  worn  by  the  pilot),  and 
compressibility  effects  of  the  flow  (such  as  shock  waves)  should  be 
assessed  as  to  their  first  or  second  order  influences  on  the  factors 
that  contribute  to  windblast  injuries. 
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SECTION  I 


INTRODUCTON 


Statistical  data  accumulated  from  five  nations  has  brought  to  light 
the  devastating  effects  of  high-speed  ejections  from  aircraft  (Glaister, 
1975).  That  is,  the  data  show  that  windblast  forces  increase  with  air- 
craft speed  to  the  point  where  an  otherwise  5-10  percent  limb-flail 
injury  rate  rises  to  40  percent  or  more.  This  is  far  from  negligible, 
but  control  of  the  forces  that  produce  excessive  motion  of  the  limbs  of 
an  ejection  seat  occupant  can  only  be  achieved  if  we  increase  our  under- 
standing of  the  aerodynamic  loading  to  which  a pilot  is  exposed  during 
high-speed  ejections.  In  a recent  report  (Schneck,  1976)  it  was  shown 
that  it  is  feasible  to  formulate  a mathematical  model  which  can  predict 
such  aerodynamic  loading.  Results  from  this  model  can  then  be  incor- 
porated into  the  Aerospace  Medical  Research  Laboratory's  modified  Calspan 
Model  of  the  Articulated  Total  Body  (ATB)  in  order  to  assess  the  kinematics 
of  limb  motion  under  the  action  of  specified  aerodynamic  forces.  In  the 
earlier  report,  potential  flow  solutions  were  presented  for  estimating 
the  pressure  distribution  around  the  forearm  of  a human  body  subjected 
to  windblast.  The  work  which  follows  expands  these  preliminary  results, 
and  extends  them  to  include  some  effects  of  flow  separation. 

The  body  model  used  here  for  the  analysis  of  windblast  forces  is 
represented  again  by  the  15-linkage  system  of  spherical,  circular- 
cylindrical,  truncated-conical  and  flat-plate  segments  shown  in  Figure  1 
(see  Schneck,  1976,  for  a complete  description).  Moreover,  specific 
attention  is  focussed  on  the  events  which  take  place  where  a body 
segment  is  either  in  contact  with  a restraining  surface,  such  as  the 
forearm  in  contact  with  an  arm  rest,  or,  equivalently,  where  two  body 
segments  are  in  contact  with  one  another,  such  as  the  upper  arm  pressing 
against  the  thorax.  For  this  situation,  the  complex  velocity  potential 
for  the  aerodynamic  cross-flow  over  the  limbs  (in  the  absence  of  flow 
separation)  has  been  shown  to  be  given  by: 

w = iraV  coth  (1) 

which  corresponds  to  the  cross-flow  streamline  pattern  depicted  in 
Figure  2 for  the  cross-sectional  geometric  configuration  shown,  with  x = 
y + iz  and  w = $ + if.  The  streamlines  illustrated  are  non-dimensionalized 
with  respect  to  V and  a. 

In  the  section  which  follows,  the  forces  tending  to  dislodge  the 
limbs  from  one  another,  or  from  a restraining  surface,  are  calculated  in 
the  absence  of  flow  separation,  i.e.,  for  the  situation  depicted  in 
Figure  2.  In  subsequent  sections,  the  results  are  then  modified  to  allow 
for  the  inclusion  of  the  more  realistic  physical  effects  of  separation 
of  flow  around  the  blunt  body  segments. 
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SECTION  II 


LIMB  DISLODGING  FORCES  IN  THE  ABSENCE  OF  FLOW  SEPARATION 


Let  p designate  the  aerodynamic  pressure  existing  at  a point  along 
the  cylindrical  surfaces  illustrated  in  Figures  2 and  3.  If  y defines 

the  azimuthal  coordinate  such  that  tan  y = (non-dimensionalized 

y 

radius  a = 1 in  Figure  2),  then  the  differential  surface  element  a(dy) 
having  unit  length  along  the  axis  of  the  limb  has  a force  on  it  of 
magnitude  df  = p(a)(dy)  directed  radially  inward  (see  Figure  3).  The 
component  of  this  force  tending  to  push  the  cylinders  together  (or  pull 


them  apart)  is  thus  given  by  dfz  = -(p  sin  y) (a  dy) , so  that  the  net 
vertical  force  becomes: 


,3ir/2 


f = / -ap  sin  y dy  per  unit  length. 

Z -tt/2 


(2) 


Now,  the  pressure,  p,  was  obtained  earlier  (Schneck,  1976)  in  the  form 
of  a pressure  coefficient  which,  for  the  cvlindrical  configuration 
illustrated  in  Figure  2,  is  given  by: 


C = 


P 1 „ 2 * 

2 pUo 


= sin  a 


2 4 / 

a 7T  , 4 ly 

7T  sech  -7T- 

2 2z 


4z 


(3) 


Before  substituting  equation  (3)  into  equation  (2)  and  performing  the 
indicated  integration,  it  is  convenient  to  make  some  transformations  and 
redefine  certain  variables.  Referring  to  Figure  3,  we  note  that  the 

2 

equation  of  the  circle  is  r = 2az,  while  r = 2a  sin  6 and  y = r cos  6. 
Thus, 


iry  ir(2a  sin  6)  (cos  9)  it  „ 

2z  = ;~2"~T~ 1 = 2 COt  6 = 6’  and’ 

4a  sin  0 


a 

IT 


(4) 


. „ ..  2 „ it  d9 

dfi  = - -j  cosec  0 d6  = - y j- 


sin  6 


2 2 2 

„ . , . z - a z , r ,4a  sin  9 , „ . 2 _ , 

Furthermore,  sin  y = 1 = — — - 1 = ~ 1 = 2sin  0-1 

3 3 2a  laL 


2 2 2 2 

sin  0 + (sin  0 - 1)  = sin  0 - cos  0 = -cos20 


so. 


- 20  + ft. 


and  dy 


2 d0. 


(5) 


Coordinate  system  defining  geometric  parameters  of  interest 


Finally, 


(6) 


And, 


2 4 2 4 

a tt  _ a tt 

4z 2 4r^ 

4a2 


2 2 
sin  0 - cos  6 


4 4 
a tt 


TT3  1 


sin2  0 


16a^sin^  9 16  sin2  0 sin2  0 


1 - cot2  0=1-^ 


8 sin2  0 


d6 

d0 


(7) 


Putting  equations  (3)  through  (7)  into  (2)  gives:  (Limits  are  on  6) 


J -(a) 


12  2 
4 pU  sin  a 
2 o 


1 " *8> 


sin2  0 


d6  ,4  . 
d0  sech  6 


+ P, 


Now 


[-cos  20 ]2  d0  (8) 

2 


, J -[y  pUQ2sin2  a + pQ] [-cos  20 ] 2 d0  = 0,  so,  with  qQ  = ^ pUQ2sin 


and  the  use  of  equation  (7),  equation  (8)  simplifies  to: 

3 , , „2  „+“> 

2(aq/ 

+<» 


/*"°°  3 , , ^2  -+°°  ao  2 

f z = J 2(aqQ^g)  sech  6 [ — j “ 1]  d<5  = -aqQTT  J sech  <5[5  - -^-J  dd  (9) 


Equation  (9)  contains  integrals  which  are  special  cases  of  the  general  form 
(see  Gradshteyn  and  Ryzhik,  page  124): 


/ 


5mdd  _ md^^cosh  6 + (n-2)dmsinh  5 
cosh11  6 (n-1)  (n-2)coshn  3 6 


n-1  J 


dmdd 


cosh11  2 6 


m(m-l)  f dm  2d6 

(n-1) (n-2)  J coshn-2 


(10) 


In  the  first  integral  of  equation  (9),  m = 2 and  n = 4.  In  the  second 
integral,  n is  again  equal  to  4,  but  this  time  m = 0.  Noting  that  the 
integrals  are  symmetric,  we  may  write,  for  the  case  m = 2: 


. 


The  first  term  on  the  right  hand  side  of  equation  (11)  evaluates  to 
zero  by  application  of  L'Hospital’s  rule.  The  second  integrates  to 

I 00 

tan  6 | = 1 (see  Gradshteyn  and  Ryzhik,  page  99)  and  the  third  integral 

°o  2k. 

is  of  the  form  (Gradshteyn  and  Ryzhik,  page  353) : / 2 

o cosh  b6 


(22k  - 2)tt21c 


b(2b) 
number 


2k 


B_.  , which,  for  k=l=b  evaluates  to  — with  the  Bernoulli 

'2k'  12 


|B„ I = — . Putting  all  of  this  together,  the  first  integral  of 
Zb  2 

equation  (9)  may  be  evaluated  in  closed  form  to  yield  (it  - 6)/9.  Now, 
for  the  case  m = 0,  equation  (10)  becomes: 


J sectJ  6 d6 
o 


sinh6  , _2  f d6 

^ ,3  Jo  3 j . 2 . 

3cosh  o o cosh  o 


(12) 


As  before,  the  first  term  on  the  right  hand  side  of  equation  (12) 
evaluates  to  zero,  and  the  second  term  integrates  to  unity.  Thus,  the 
second  integral  of  equation  (9)  may  be  evaluated  in  closed  form  to  yield 
2 

(4/3)  (-it  /4).  We  may  then  write: 

2 , .2  -aq  ir  „ 

fz  = -av[JLir  - = -t~  ["2lt  - 6]  = 9aq0  (13) 

Note  the  positive  resultant  sign  of  equation  (13) , indicating  that  the 
net  force  per  unit  length  of  limb  is  in  the  positive  z-direction,  or 
tending  to  cause  dislodgement  of  the  limb  from  a restraining  surface  (which 
could  also  be  another  limb) . 


In  order  to  evaluate  equation  (13)  further,  we  need  to  say  something 
about  q and  a.  As  a reasonable  approximation,  we  may  begin  by  assuming 
that  the  radius  of  the  human  forearm  is  on  the  order  of  1.5  inches  (average) 
and  that  the  length  of  this  portion  of  the  arm  is  about  one  foot.  Thus, 
the  total  force  acting  on  the  forearm  can  be  approximated  by  the  relation 
9 

f - -«■  q . Now,  consider  an  ejection  taking  place  at  an  altitude  of 

t o O 

10,000  feet.  At  this  height,  the  density  of  air  relative  to  a known 


standard  is  given  by  (Shames,  page  540)  = 0.7385,  where  Po  = 0.002378 

3 ° 

slug/ft  , and  the  speed  of  sound,  c = 1,078  feet  per  s cond.  Thus, 

lp  2222  22 

q = — — PQ[U0  /c  ]c  sin  a = 1020M  sin  a,  where  M < >rresponds  to  the 

Ko 

Mach  Number  of  the  flow.  Substituting  this  result  into  the  equation  for 
f gives  the  total  limb-dislodging  force  as  a function  of  Mach  Number 
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and  Angle  of  Attack  at  the  specified  altitude  of  10,000  feet.  The  final 
equation  is: 

f = 1148  M2  sin2a  (14) 

Figure  4 shows  equation  (14)  plotted  as  limb  dislodging  force  vs. 

Mach  Number  with  the  angle  of  attack  appearing  as  a parameter.  Although 
the  figure  illustrates  force  behavior  in  che  supersonic  range  as  well  as 
the  subsonic,  it  is  recognized  that  the  results  past  M = 1 are  not  realistic, 
since  the  analysis  has  thus  far  neglected  the  presence  of  any  shock  waves 
in  the  flow.  In  any  case,  the  subsonic  values  for  f are,  in  themselves, 
very  revealing.  For  example,  at  the  higher  angles  or  attack  one  observes 
a rapidly  increasing  limb-dislodging  force.  The  significance  of  these 
large  forces  may  be  appreciated  by  examining  them  in  relation  to  the 
super-imposed  ordinate  axis  which  is  a measure  of  several  published  pro- 
bability results.  That  is,  in  a study  conducted  by  Horner  and  Hawker  (1973) 
it  was  found  that  a pilot's  average  grip  retention  is  such  that  the  pro- 
bability of  his  letting  go  is  100%  if  the  dislodging  force  exceeds  some 
600  pounds.  From  Figure  4,  we  see  that  f = 600  pounds  when  M = .72  for 
a = 90°  and  when  M = .83  for  a = 60°.  Thus,  at  these  higher  angles  of 
attack  a pilot's  musculo-skeletal  system  is  not  likely  to  withstand  the 
tendency  for  dislodgement  from  a restraining  surface  if  he  is  ejecting 
at  Mach  Numbers  in  excess  of  around  0.7.  Similarly,  in  a study  conducted 
by  Payne  (1975)  it  was  found  that  the  probability  of  major  flail  injury 
is  around  100%  if  the  ejection  Mach  Number  exceeds  1.25.  This  would 
appear  to  be  self-consistent  with  the  predicted  limb  dislodging  forces 
as  interpreted  above  together  with  the  results  of  Horner  and  Hawker.  To 
this  level  of  approximation,  therefore,  one  would  have  to  conclude  that 
the  generation  of  stagnation  points  in  the  flow  produces  forces  that  can 
cause  limb-dislodgement  (with  subsequent  flail  and  possible  serious 
injury).  Moreover,  these  forces  are  sensitive  to  the  angle  at  which  the 
limb  intercepts  the  flow,  such  that  the  higher  the  angle,  the  greater 
the  tendency  for  dislodgement.  And  finally,  the  forces  increase  rapidly 
with  speed  of  ejection,  which  correlates  well  with  the  finding  that 
windblast  injuries  increase  dramatically  as  a function  of  airspeed 
(Glaister,  1975).  Figure  5,  taken  from  Payne,  Hawker  and  Euler  (1975) 
is  a photograph  of  an  ejection  seat  occupant  sitting  in  the  ACES-II  Seat 
at  -15°  Yaw  and  -15°  Pitch,  during  a wind-tunnel  simulation  of  an 
ejection  from  the  F-105.  The  arrows  point  out  several  critical  regions 
of  the  pilot-ejection-seat  configuration  where  flow  stagnation,  with  its 
consequent  limb-dislodging  force  distribution,  is  likely  to  occur. 

Observe,  in  particular,  those  regions  where  the  pilot  is  gripping  the 
seat,  where  his  upper  arms  are  in  contact  with  his  torso,  where  the 
lower  legs  are  in  contact  with  the  seat  pan  and  where  his  head  is  in 
contact  with  the  back  of  the  seat.  All  of  these  represent  sites  of 
potentially  serious  windblast  and  flail  injury.  Based  on  these  results, 
one  is  therefore  encouraged  to  pursue  the  mathematical  analysis  with  a 
certain  degree  of  assurance  that  the  approach  is  faithfully  describing 
events  as  they  have  been  experimentally  observed  to  take  place.  With 
this  in  mind,  we  now  proceed  to  modify  the  theory  to  include  flow 
separation  effects. 
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VORTEX  MOTION  IN  A REGION  OF  SEPARATION 


In  aerodynamic  flow  theory,  two  simplifications  have  been  conveniently 
and  commonly  employed  to  investigate  separated  flows  over  arbitrary 
three  dimensional  bodies  of  revolution  (Marshall  and  Deffenbaugh,  1975). 

The  first  asserts  that  there  is  a direct  analogy  between  three-dimensional 
steady  flow  and  two-dimensional  unsteady  flow  (Allen  and  Perkins,  1951)  — 
such  that  a three-dimensional  steady  separated  flow  problem  can  be 
analyzed  as  a two-dimensional  unsteady  separated  flow  problem.  The 
second  is  based  upon  the  observation  (verified  experimentally)  that  a 
two  dimensional  unsteady  wake  can  be  described  by  a distribution  of 
inviscid  vortices,  originating  from  the  separation  of  shear  layers  and 
superimposed  on  the  unseparated  potential  flow  solution  (Mello,  1959, 
Sarpkaya,  1968  and  Marshall  and  Deffenbaugh,  1975).  These  vortices  are 
modified  by  diffusion.  It  is  the  second  of  these  two  assumptions  that  we 
shall  exploit  here,  since  arguments  have  already  been  presented  for 
considering  flow  over  the  forearm  in  contact  with  a restraining  surface 
to  be  two-dimensional  (Schneck,  1976) . 

Thus,  to  examine  the  separated  flow  patterns  around  the  geometric 
configuration  illustrated  in  Figure  2,  consider  the  situation  depicted 
in  Figure  6.  A Vortex  pair  has  been  placed  on  the  lee  or  downstream 
side  of  the  two  cylinders  as  shown.  Neither  the  location,  yQ  + izQ,  nor 

the  strength,  k,  of  these  two  vortices  can  be  specified  at  this  time. 

In  fact,  the  location  of  the  vortex  pair  can  not  be  determined  uniquely 
from  potential  flow  theory  alone  (see  later  discussion).  However,  the 
vortex  strength  can  be  computed  if  we  say  something  about  how  these 
vortices  move  relative  to  the  cylinders.  In  this  respect,  reasoning 
developed  in  the  earlier  report  (Schneck,  1976)  has  suggested  that, 
during  the  first  few  milliseconds  of  high-speed  flow  development  around 
a blunt  body  of  revolution,  the  separated  wake  stays  confined  to  the 
immediate  region  of  the  body,  and  consists  essentially  of  two  thin 
vortex- layers,  symmetrically  situated.  One  may  thus  model  this  separated 
flow  analytically  by  super- imposing  stationary  vortex  filaments  upon  the 
linear  inviscid  flow  solutions  obtained  earlier.  The  condition  of 
stationarity  then  allows  the  vortex  strength  to  be  determined  as  described 
later  in  this  report. 

Derivation  of  Complex  Velocity  Potential  for  a Vortex  Pair 
Downstream  of  Two  Tangent  Cylinders  with  Parallel  Axes 

In  order  to  super- impose  the  vortex  layers  on  the  earlier  solutions, 
as  described  above,  we  must  first  develop  a complex  velocity  potential 
to  describe  the  vortex  pair  illustrated  in  Figure  6.  The  analytic  form 
for  this  potential  can  be  deduced  in  part  by  examining  the  transformation 
that  maps  the  double  cylinder  configuration  of  Figure  6 into  a single 
unit  circle  centered  at  the  origin  of  the  transformed  plane  (see  Carrier, 
Krook  and  Pearson,  1966,  pgs.  132-134).  Consider  the  bilinear  trans- 
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formation: 


u + i 7 


*X  ~ 2 , ixx  2 2X 

X XX 


i(y2  + z2)  - 2(y  - iz) 

2 2 
y + z 


(15) 


On  the  upper  cylinder  of  Figure  6,  with  a ■ 1,  the  equation  of  the 
2 2 

circle  is  y + z = 2z.  Substituting  this  into  equation  (15),  we  find 
that  v = 2 » constant  for  any  value  of  u.  Similarly,  on  the  lower 
2 2 

cylinder,  y + z = -2z,  yielding  v = 0 for  any  value  of  u.  The  bilinear 
transformation  (15)  thus  maps  the  two  tangent  circles,  S^  and  S2  of 

Figure  6 into  the  parallel  lines  v » 2 and  v = 0,  respectively,  in  the 
u,  v plane.  This  is  illustrated  in  Figure  7A. 

Next,  consider  the  exponential  transformation: 

J (u  + iv)  TU  Tiv  J u 

x + is  « e =e  e = e (cos  j v + isin  j v) . (16) 


For  v - 2,  s » 0 and  x » which  maps  S^  onto  the  negative  real 

axis  of  the  x,s  plane.  Likewise,  for  v * 0,  s = 0 and  x = 

which  maps  S^  onto  the  positive  real  axis  of  the  x,s  plane.  The 

exponential  transformation  (16)  thus  maps  the  two  parallel  lines  of 
Figure  7A  onto  the  real  axis  of  the  x,s  plane,  as  illustrated  in  Figure 


Finally,  consider  the  inverse  bilinear  transformation: 

x + is  + i x + i (x  + i)  (x  + i)  _ x2  + 2ix  - 1 
4 n " X + is  - i " x - i " (x  - i)(x  + i)  = x2  + 2 * 


(17) 


where  s has  been  set  equal  to  zero  in  accordance  with  the  discussion  of 

2 2 

equation  (16).  From  equation  (17)  we  get  £ = (x  - l)/(x  +1)  and 


4 2 2 

„ ,,  2 . „ ,2  . 2 x - 2x  + 1 + 4x 

n - 2x/(x  + 1).  Thus,  £ + n * 5 5 

(x  + 1) 


x*  + 2x2  + 1 
x4  + 2x2  + 1 


1, 


so  that,  for  s = 0 equation  (17)  defines  the  equation  of  a circle  of  radius 
1 in  the  £,n  plane.  The  inverse  transformation  (17)  thus  maps  the  real  axis 
of  the  x,s  plane  into  a unit  circle  centered  at  the  origin  in  the  £,n  plane, 
as  illustrated  in  Figure  7C. 
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Now,  putting  together  the  transformations  (15),  (16)  and  (17),  one 
obtains  the  single  function  that  maps  the  two  unit  cylinders  of  Figure  7A 
into  the  single  unit  circle  of  Figure  7C.  Thus, 


W = 5 + ih 


x + is  + i 
x + is  - i 


§ (u  + iv) 
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\ (u  + iv) 
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TT  . 
21 


e 


(18) 


By  multiplying  the  top  and  bottom  of  the  quotient  in  equation  (18)  by 
+ (tt  / 2y) 

the  quantity  e , the  mapping  function  may  be  put  into  the  form: 


W = 


e 


TT 

2X 


+ 

+ e 


TT 

2x 


- coth 


(— ) 
C2X;‘ 


(19) 


Figure  8 illustrates  the  net  transformation  represented  by  equation  (19). 
The  significance  of  this  transformation  is  that  it  suggests  an  analytic  form 
for  the  complex  potential  when  inviscid  flows  encounter  boundaries  such  as 
those  depicted  in  Figures  2,  6 or  8.  For  example,  let  us  expand  equation 
(19)  into  its  real  and  imaginary  parts  and  see  what  information  can  be 
gained  from  the  results  obtained.  Expansion  of  the  hyperbolic  cotangent 
of  complex  argument  yields  the  following: 


coth  (^) 


cothfa^  - coth  "fr  - “>■  - 
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r r 


A1(y,z)  + iA2(y,z)  = A3(y,z)elg^y’Z^ 


(20) 


sinh 


where:  A^(y,z) 
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/~2  2 A (y,z) 

A (y,z)  = /A  (y,z)  + A„  (y,z);  and,  B(y,z)  = Arc  Tan  7—7 7 

J 1 1 A1(y,z) 

2 

Now  observe  that  if  r = 2z  (the  equation  of  the  upper  cylinder  in  Figures 
2 or  8)  then  A^y.z)  = (sinh  ~^)/(cosh  and  A2(y,z)  = (-l)/(cosh  , so 


2 o o (sinh  + 1 cosh2  ^ 

that:  A„  (y,z)  = A..  (y,z)  + A„  (y,z)  = - r — = - — - = 1. 

-3  J.  2 / , Try ^ 2 , 2 Try 

(cosh  -r^)  cosh  -t~- 

2z  2z 


The  exact  same  statement  can  be  made  for  the  lower  cylinder  which  is  defined 
2 

by  the  equation  r = -2z  because  the  negative  sign  disappears  when 
A^(y,z)  and  A,,(y,z)  are  squared  in  the  calculation  of  A^ty.z).  The  modulus 

of  equation  (20)  is  thus  equal  to  1 on  the  boundary  defined  by  figures  2 
or  8.  Moreover,  the  logarithm  of  this  modulus  is  therefore  equal  to  zero 
on  this  particular  boundary.  Consequently,  we  conclude  that  the  quantity 
represented  by  the  logarithm  of  the  modulus  of  equation  (20)  could  represent 
the  stream  function  for  flow  around  a double-cylinder  geometric  configuration 
such  as  that  being  considered  here.  In  fact,  upon  plotting  this  logarithm, 
we  find  that  it  corresponds  to  a circulation  pattern  around  just  such  a 
double  cylinder  configuration.  This  is  illustrated  in  Figure  9. 

To  summarize,  up  to  this  point  it  has  been  shown  that  if  we  take  the 
log  of  both  sides  of  equation  (20),  the  real  part  of  the  right  hand  side, 
i.e.,  log  A^(y,z)  represents  the  streamline  pattern  for  fluid  circulation 

around  a boundary  consisting  of  two  unit  cylinders  tangent  to  one  another 
as  shown  in  Figures  2,  8 and  9.  This  log  function  should  therefore 
correspond  to  the  stream  function,  of  the  complex  velocity  potential 
suggested  by  equations  (19)  and  (20) . In  order  to  affect  this  corres- 
pondence we  must  multiply  the  log  of  equation  (20)  by  the  quantity  i,  so 
that  log  A-j(y,z)  will  now  become  the  imaginary  part  of  the  complex  potential 

w = 4>  + iV,  and  thus  correctly  represent  a stream  function.  Finally,  to 
designate  some  strength  to  the  circulation  pattern,  we  assign  some  arbitrary 
parameter  K to  the  complex  velocity  potential  giving  the  result: 


w = iK  log  [coth  (-j^)]  (21) 

The  imaginary  part  of  equation  (21)  now  yields  a circulation  streamline 
pattern  (Figure  9)  around  the  boundaries  defined  by  the  equation 

2 

r = + 2z.  This  circulation  pattern  is  a measure  of  the  strength  of  vortices 
which  are  downstream  of  the  boundaries,  and,  in  fact,  is  directly  propor- 
tional by  a factor  of  2tt  to  the  magnitude,  k,  or  vortex  strength  of  the 
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Figure  9.  Circulation  Pattern  around  a Double— Cylinder  Configuration. 

The  Non-dimensionalized  Stream  Functions  for  this  Circulation 
Pattern  are  given  by  the  logarithm  of  Equation  (20)  defined 
in  the  text. 


vortex  pair  defined  by  the  complex  velocity  potential: 


Thus,  we  conclude  that  the  vortex  pair  on  the  downstream  side  of  the 
cylinders  illustrated  in  Figure  6 may  be  defined  by  the  stream  function: 


x — x 

Y = ^ Im  UK  log  [coth  )]}  Im  [i  log  --]  (23) 

x X - xD 

Note,  in  equation  (23),  that  K (a  real  constant)  and  the  location,  y + iz 

o — o 

are  still  undefined.  In  any  event,  we  may  plot  (2ttT/K)  for  arbitrary  y and 

ZQ  just  to  get  an  idea  of  what  the  streamline  patterns  look  like.  This  is 

done  in  Figure  10  for  a vortex  pair  placed  at  y =1.6  and  z = + 1.6. 

o o — 

Determination  of  Circulation  Strength 

It  was  mentioned  earlier  that  the  strength,  K,  of  the  circulation 
pattern  could  be  computed  if  we  say  something  about  how  the  vortex  pair 
moves  relative  to  the  cylinders.  More  specifically,  reasoning  has  been 
presented  for  assuming  that  early  in  the  flow-development  process  the 
separated  wake  stays  confined  to  the  immediate  region  of  the  cylinders, 
and  that  the  vortex  layers  therefore  remain  stationary  with  respect  to 
these  cylinderical  bodies.  Stating  this  fact  mathematically  allows  K to 
be  determined  as  follows:  suppose  the  flow  we  are  examining  consists  of 
a cross-flow  pattern  having  streamlines  defined  by  the  imaginary  part  of 
equation  (1),  and  vortices  defined  by  the  stream  functions  of  equation 
(23).  Then,  the  cross-flow  would,  for  K arbitrary,  have  a tendency  to 
carry  the  vortex  pair  along  with  it,  so  that  these  would  be  shed  down- 
stream as  a function  of  time.  Moreover,  each  vortex  would,  itself,  have 
a tendency  to  induce  a translational  velocity  at  the  center  of  its 
corresponding  partner.  In  fact,  for  the  situation  shown  in  Figure  6, 
note  that  the  velocity  field  induced  by  each  vortex  on  the  other  has  a 
tendency  to  cause  the  vortex  pair  to  drift  upstream.  Now,  if  we  impose 
the  stationarity  condition  as  described  earlier,  then,  by  requiring  the 
downstream  velocity  of  the  vortex  pair  to  exactly  balance  the  upstream 
velocity,  we  find  that  K no  longer  becomes  arbitrary  but  takes  on  a 
value  determined  uniquely  by  the  condition  of  stationarity.  The  velocity 
induced  at  y^  + iz^  by  the  complex  velocity  potential  given  in  equation 

(1)  is  (dw/dx)  evaluated  at  y , zq.  Similarly,  the  velocity  induced 

at  yQ  + izQ  by  the  stream  function  given  in  equation  (23)  is 
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(— 3^ / 9 z ) + iO’f/Sy)  also  evaluated  at  y , zq.  Adding  these  two  velocities 

together  and  setting  their  sum  equal  to  zero  thus  insures  stationarity 
of  the  vortex  pair  and  gives  a relationship  between  the  circulation 
strength,  K,  the  cross-flow  velocity,  V,  and  the  vortex  locations  y^  + 

izQ.  At  this  point,  however,  we  cannot  pursue  the  analysis  any  further 

until  yQ  and  zq  are  somehow  quantitatively  defined. 

Potential  flow  theory  alone  can  provide  no  information  concerning 
the  location  of  vortices  in  a separated  wake.  Such  information  must  be 
determined  either  by  considering  the  variation  of  pressure  across  the 
boundary  layer,  as  described  by  Sarpkaya  (1968)  and  Bar-Lev  and  Yang 
(1975),  or  by  performing  a series  of  carefully  designed  experiments.  It 
is  our  intent  to  exploit  both  of  these  techniques  in  future  studies. 


SECTION  IV 
CONCLUDING  REMARKS 


In  this  report,  the  mathematical  model  developed  in  an  earlier  work 
(Schneck,  1976)  has  been  expanded  upon  in  two  ways.  First,  actual 
forces  which  tend  to  dislodge  two  limbs  from  one  another,  or  from  a 
restraining  surface,  have  been  calculated  in  the  absence  of  flow  separation. 
Second,  the  theory  has  been  modified  to  allow  for  the  inclusion  of  some 
physical  effects  attributable  to  separation  of  flow  around  the  blunt 
body  segments.  It  is  clear  from  the  results  obtained  thus  far,  that 
aircraft  ejections  taking  place  above  Mach  0.8  have  associated  with  them 
a very  high  probability  of  limb-flail  injury.  This  would  appear  to  be 
due  to  the  forces  generated  as  a consequence  of  the  manner  in  which  the 
ejection  seat  occupant  intercepts  the  air  stream.  The  generation  of 
stagnation  points  in  the  flow  leads  to  limb-dislodging  forces  which, 
beyond  Mach  0.7,  exceed  the  musculo-skeletal  ability  of  the  pilot  to 
"hold  on".  This,  then,  becomes  an  important  design  criterion  for 
ejection  systems,  i.e.,  to  minimize  stagnation  of  the  flow  as  the  pilot 
intercepts  the  air  stream. 

It  is  likely  that  flow  separation  also  plays  an  important  role  in 
the  generation  of  limb  dislodging  forces  because  this  phenomenon  is 
responsible  for  the  major  component  of  drag  when  a blunt  body  is  placed 
in  a high-speed  fluid  flow.  It  is  thus  important  to  quantify  the 
effects  of  flow  separation  and  it  is  anticipated  that  this  will  be  done 
in  future  investigations  — both  analytically  and  experimentally.  Also 
to  be  examined  in  future  studies  is  the  time-course  of  limb  dislodging 
forces.  That  is,  the  force  distribution  shown  in  Figure  4 exists  at  the 
onset  of  ejection  and  these  are  peak  forces  that  prevail  so  long  as  the 
limb  is  in  contact  with  a restraining  surface.  Once  contact  with  the 
surface  is  broken,  there  follows  a redistribution  of  forces,  and  this 
information  is  desirable  input  to  the  ATB  simulation  which  provides 
kinematic  data.  All  of  this  work  is  being  carried  out  in  an  effort  to 
define  specific  design  criteria  for  safer  ejections. 


LIST  OF  SYMBOLS 


A^y.z) 


A2(y,z) 


A3(y,z) 


C 

P 


F 

K 


M 


R 


o 


sinh  ^ 
z 

r 

. Try  trz 

cosh  2 “ cos  — 2 
r r 

, TTZ 
-sin  — j 

r 

, Try  ttz 

cosh  — j - cos  — ^ 

r r 


v412(y,z)  + A22(y,z) 


Bernoulli  Number  of  order  2k 


Surface  Pressure  Coefficient  = 


The  same  as  f 

Circulation  Strength  around  double-circular-cylinder  configuration 


U 

Mach  Number  = — 
c 


Designates  curves  in  space 

Magnitude  of  Free-Stream  Air  Flow  Velocity 

Magnitude  of  the  Cross-Flow  Velocity 
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m 


Mapping  Function 


Radius  of  circles  tangent  to  y-axis;  Average  human  forearm  radius 


Arbitrary  Coefficient 


Speed  of  Sound 


Ordinary  Differential  Operator 


Base  for  Natural  Logarithms  = 2.718281828 


Magnitude  of  the  radial  force  generated  by  the  fluid  pressure,  p 


Total  vertical  force  tending  to  dislodge  the  human  forearm  from 
a surface  (or  other  limb)  with  which  it  is  in  contact. 


The  component  of  radial  force  in  the  z-direction,  per  unit  length 
perpendicular  to  the  y-z  plane. 


arbitrary  integer  exponent  or  subscript 


arbitrary  integer  exponent  or  subscript,  independent  of  k 


arbitrary  integer  exponent  or  subscript,  independent  of  k and  m 


used  as  a subscript  to  designate  free-stream  conditions,  or  to 
locate  specific  points  in  space,  or  to  specify  reference 
quantities. 


Aerodynamic  Pressure 


Reference  Pressure 
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12  2 
-=■  pU  sin  a 
l o 


2 2 2 

Radial  coordinate,  r = y + a 


Radial  location  of  the  center  of  a vortex 


Exponential  Transformation  Coordinate  (See  x below) 


Time  Coordinate 


Bilinear  Transformation  Coordinate  (See  v below) 


Bilinear  Transformation  Coordinate  for  the  transformation  u + iv 


Complex  Velocity  Potential,  w = <J>  + if 


Exponential  Transformation  Coordinate  for  the  transformation 
x + is 


Coordinate  axis  perpendicular  to  the  line  along  which  two  limbs 
are  in  tangent  contact,  or  the  line  along  which  one  limb  is  in 
tangent  contact  with  a restraining  surface,  and  lying  within 
the  restraining  surface. 


the  y-coordinate  locating  the  center  of  a vortex. 


Coordinate  axis  orthogonal  to  y and  the  restraining  surface, 
i.e.,  normal  to  the  arm  rest. 


The  z-coordinate  locating  the  center  of  a vortex 


The  angle  of  attack  of  the  free  stream  relative  to  the  center- 
line  of  the  limb. 


B(y,z)  Arc  Tan 


A2(y,z) 

Ax(y,zf 


M|=l 


n 

K 

P 

P, 

Y 

0 

0 

c 

a 

T 

c 


X 

X 

n 


z - a 

Azimuthal  Location  of  Vortex  = Arc  Tan  — 


Transformation  Coordinate  (See  £ below) 
Vortex  Strength 
Fluid  Mass  Density 
Reference  Density 

Azimuthal  Coordinate,  Tan  y = — ~--a 
Inscribed  Angle 

z 

Vortex  Location  = Arc  Tan  — 


z - a 


y 


o 


y + io  = re 


Inverse  Bilinear  Transformation  Coordinate  for  the  transformation 
£ + in. 


cot  0 

, . „ i0 

y + iz  = Re 

y - iz 

3.141592654 


00  Designates  Infinity 

f Stream  Function  of  the  Flow 
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Velocity  Potential  of  the  Flow 


Designates  Partial  Differentiation 


APPENDIX 


- ySO’CTT' 


COMPUTER  PROGRAMS  FOR  CALCULATING 
AND  PLOTTING  STREAM  FUNCTIONS 


A.  NON-DIMENS 10NALIZED  STREAMING  FLOW  PATTERN 


DIMENSION  G(4200) ,H(4200) ,LABEL(10) ,LABLE(10) ,N2(400) 
DATA  LABEL (1) / ' Y '/ 

DATA  LABLE (1) / ' Z '/ 

Y=.5 
X=5 . 

CALL  SAXIS (X,Y ,LABLE, -4, 7.0,90.,1.,1.,1.,1.) 

X=0. 

CALL  SAXIS (X,Y, LABEL, -4, 10. ,0. ,-5. ,1. ,1. ,0.0) 

CALL  CIRCLE(5. ,1.5,1.) 

PI=355./113. 

DO  9 N-1,31,2 
PSIB=0. 1*N 
PSIU=PSIB+. 005 
PSIL=PSIB-. 005 
N1=0 

DO  4 J=l,400 
Y= . 025* (J-l)-5 . 0 
N2(J)=0 
11=166 

IF(Y.GT.-2.0.AND.Y.LT. 2. 0)11=500 
DO  5 1=1,11 

IF(Y.GT.-2.0.AND.Y.LT. 2.0) GO  TO  7 
Z=. 030*1 
GO  TO  10 
7 Z=. 010*1 

10  IF(((Z-l)*(Z-l)+Y*Y).LT.l)GO  TO  5 
RSQ=Y*Y+Z*Z 
IF (RSQ. EQ. 0. 0)GO  TO  5 
Q=RSQ-Z/2 

IF(Q.GT.-.01.AND.Q.LT. .01) GO  TO  5 

F1Y=PI*Y/RSQ 

F1Z=PI*Z/RSQ 

IF (ABS (F1Z-1. 5708) . LT. 0. 001)GO  TO  5 
F2=TAN (F1Z) 

F3=C0SH(F1Y) 

F4=TANH(F1Y) 

PSI=PI*F2/((F3*F3)*(F4*F4+F2*F2)) 

IF (PSI . LE. 0. )G0  TO  5 
IF  (PSI.LT.PSIL)GO  TO  5 
IF  (PSI . LT. PSIU. AND. PSI . GT. PSIL)GO  TO  45 
GO  TO  5 
45  N1=N1+1 

N2(J)=N2(J)+1 
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IF(N2(J) .GE.2)  GO  TO  92 
G(Nl)=Y+5. 

H(N1)=Z+. 5 
GO  TO  5 
92  N1=N1-1 
5 CONTINUE 
4 CONTINUE 

CALL  LINE(G,H,N1,+1) 

9 CONTINUE 
CALL  PLOT(0. ,0. ,-4) 

STOP 


B. 


NON-DIMENS IONALIZED  VORTEX  PATTERN 


DIMENSION  G(4200) ,H(4200) , LABEL (10) , TABLE (10) ,GG(3000) ,HH(3000) 
1,N2(400) ,Z1(600) ,W(402) ,Z2(1500) 

DATA  LABEL(l) / ' Y '/ 

DATA  LABLE (1) / ' Z '/ 

DATA  YO/1. 6/ ,ZO/l. 6/ 

Y=.  5 
X=5. 

CALL  SAXIS(X,Y,LABLE,-.4, 7.0,90.  ,1-  .1.  ,1-  ,1.) 

X=0. 

CALL  SAXIS(X,Y,LABEL,-4, 10. ,0. ,-5. ,1. ,1. ,0.0) 

CALL  CIRCLE(5. ,1.5,1.) 

PI=355./113. 

DO  9 N-1,31,2 

PSIB=0.1*N 

PSIU=PSIB+. 005 

PSIL=PSIB-. 005 

KK=1 

N1=0 

N3=0 

DO  4 J=1 , 400 
Y=. 025* (J-l)-5. 0 
N2(J)=0 
11=166 

IF(Y. GT. 0.00. AND. Y.LT. 3. 0)11=500 
DO  5 1=1,11 

IF ( Y . GT . 0 . 00 . AND . Y . LT . 3 . 0) GO  TO  7 
Z=. 030*1 
GO  TO  10 
7 Z=. 010*1 
10  W(J)=Z 

A4A=SQRT ( (Y-YO)**2+(Z-ZO)**2) 

B4=SQRT ( (Y-YO) **2+(Z+ZO)**2) 

IF (A4A. EQ. 0. 0)PSI=1000 
IF (A4A . EQ . 0 . 0) GO  TO  5 
E4=B4/A4A 

IF(E4.EQ.0.0)PSI=1000 
IF(E4.EQ.O.O)GO  TO  5 
PSI=ALOG(E4) 

IF(PSI.LE.O.)GO  TO  5 
IF  (PSI.LT.PSIL)GO  TO  5 
IF  (PSI.LT.PSIU.AND.PSI.GT.PSIL)GO  TO  45 
GO  TO  5 
45  N1=N1+1 
Z1(1)=W(J) 

N2(J)=N2(J)+1 

IF(N2(J) ,GE.2)GO  TO  92 

IF(N. EQ. 1. AND.KK. EQ. 2) GO  TO  12 
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j . 


wmP 


r 


IF (ABS  (Z1  (KK)-W(J)  ) . C>T . 0. 30)G0  TO  92 
12  KK=KK+1 
Z1(KK)=W(J) 

G(Nl)=Y+5. 

H(N1)=Z+. 5 
GO  TO  5 
92  N1=N1-1 

IF (ABS (Z 1 (KK) -W ( J ) ) . LE . 0 . 05 ) GO  TO  5 

N3=N3+1 

Z2(l)=1000 

IF (ABS (Z2 (N3)-W(J) ) . LE. 0. 05)GO  TO  8 
Z2(N3)=W(J) 

GG(N3)=Y+5 
HH(N3)=Z+. 5 
GO  TO  5 

8 N3=N3-1 
5 CONTINUE 
4 CONTINUE 

CALL  LINE(G,H,N1,+1) 

IF(N3.LE.2)GO  TO  9 
CALL  LINE(GG,HH,N3,+1) 
WRITE(6,50)N1,N3,PSIB 
50  FORMAT (6X, 14, 2X, 14, 2X,F10. 3) 

9 CONTINUE 

CALL  PLOT (0. , 0. ,-4) 

STOP 

END 


34 


C.  NON-DIMENS IONALIZED  CIRCULATION  PATTERN 


DIMENSION  G(4200) ,H(4200) , LABEL (10) , LABLE (10) ,N2(500) 
DATA  LABLE (1) / ' '/ 

DATA  LABEL(l) / ' '/ 

Y=0. 5 
X=5. 

CALL  SAXIS  (X,Y, LABLE, -4, 10. ,90. ,-5. ,1. ,1. ,0. ) 

X=1.0 

Y=5. 5 

CALL  SAXIS (X, Y,LABEL,-4, 8.0, 0.0, -4. ,1. ,1. ,0.) 

CALL  CIRCLE(5. ,4.5,1.) 

CALL  CIRCLE(5. ,6.5,1.) 

PI=355./113. 

DO  15  L=l,2 
F=l.  0 

IF (L.EQ. 2)F=-1. 0 
DO  9 N=7 , 152 , 15 
PSIB=N* . 005 
PSIL=PSIB-.01 
PSIU=PSIB+. 01 
N1=0 

DO  4 J=l, 350 
Y=0. 02*(J-l)-3. 5 
N2(J)=0 
11=120 

IF (Y. GT. -1.2. AND. Y.LT. 1.2) 11=300 
DO  5 1=1,11 

IF (Y.GT. -1.2. AND. Y.LT. 1. 2) GO  TO  7 
Z=. 025*1-. 025 
GO  TO  10 
7 Z=. 010*1-. 010 
10  Z=F*Z 

IF (Y . EQ . 0 . 0 . AND . Z . EQ . 0 . 0) GO  TO  5 
IF(((Y+1)*(Y+1)+Z*Z) .LT.l)GO  TO  5 
IF(((Y-l)*(Y-l)+Z*Z).LT.l)GO  TO  5 
R=SQRT(Y*Y+Z*Z) 

DP=PI*Z/R/R 

IF(ABS(DP).GT.150)GO  TO  5 
Al=COSH (PI*Z /R/R) -COS (PI*Y /R/R) 

IF (ABS (Al) . LT. 0. 0001)GO  TO  5 
A2=SINH(PI*Z/R/R) /Al 
IF(ABS(A2) .LT. 0.001) GO  TO  5 
A3=0. 0-SIN (PI*Y/R/R) /Al 
IF (ABS (A3) .LT. 0.001) GO  TO  5 
A4=SQRT (A2*A2+A3*A3) 

IF (A4. EC. 0.0) GO  TO  5 
PSI=ALOG(A4) 

IF(PSI.LE.O.)GO  TO  5 
IF  (PSI.LT.PSIL)GO  TO  5 
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IF  (PSI.LT.PSIU.AND.PSI.GT.PSIL)GO  TO  45 
GO  TO  5 
45  N1=N1+1 

N2  (J)=N2 (J)+l 
IF(N2(J) ,GE.2)GO  TO  2 
G(Nl)=z+5 . 0 
H(Nl)=y+5. 5 
GO  TO  5 
2 N1=N1-1 
5 CONTINUE 
4 CONTINUE 
CALL  LINE(G,H,N1,+1) 

9 CONTINUE 
15  CONTINUE 

CALL  PLOT(0. ,0. ,-4) 

STOP 

END 
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